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Abstract 


This  paper  treats  the  problem  of  estimating  the  intensity 
distribution  of  an  optical  pattern  which  has  been  distorted  by  diffraction 

1 

i 

3 

and  noise.  Under  the  assumption  that  the  signal  and  noise  are  additive 
and  that  they  have  prescribed  means  and  covariance  matrices,  the 

3 

optimum  linear  estimate  of  the  object  distribution  is  obtained.  This 
estimate  is  shown  to  be  a  linearly -filtered  combination  of  the  a  priori 
mean  and  the  observed  image  distribution,  with  the  details  of  the  filtering 
depending  on  the  prior  information,  the  noise  statistics,  and  the  optical 
imaging  system.  The  performance  of  the  optimum  estimation  procedure 
is  evaluated  for  the  important  special  case  of  large  a  priori  uncertainty. 

Finally,  the  optimum  procedure  for  processing  multiple  observations 
is  discussed,  ft  is  shown  that  as  the  number  of  observations  becomes 
large,  the  estimate  approaches  the  true  object. 


4r 


I,  Introduction 

The  problem  of  restoring  the  detail  removed  from  an  optical  pattern 
by  diffraction  has  been  considered  by  Harris  [  l],  Using  some  well- 
known  results  from  the  theory  of  Fourier  transforms,  he  showed  that 
in  the  absence  of  noise,  an  optical  pattern  known  a  priori  to  be  of  bounded 
spatial  extent  can  be  exactly  reconstructed  despite  the  rejection  of  high 
spatial  frequencies  by  an  optical  system. 

While  the  noise-free  problem  is  of  considerable  theoretical  interest, 
it  does  not  accurately  represent  the  physical  world.  In  any  actual 
situation,  spurious  disturbances  such  as  background  radiation,  film 
granularity,  and  detector  noise  are  superimpose!  on  the  quantities  to  be 
measured.  Thus,  any  restoration  procedure  will  involve  a  certain 
amount  of  unavoidable  error.  Furthermore,  a  given  restoration  procedure 
may  perform  very  badly  in  the  presence  of  noise  even  though  it  repro¬ 
duces  the  pattern  exactly  when  noise  is  absent.  Hence,  it  is  important 
to  take  into  account  any  knowledge  we  have  ab^ut  the  noise  statistics,  as 
well  as  the  a  priori  information  about  the  pattern  to  be  estimated,  in 
designing  a  restoration  procedure. 

In  this  paper,  we  treat  the  problem  of  obtaining  optimum  estimates 
of  an  optical  object  in  the  presence  of  additive  noise  using  the  criterion 
of  minimum  mean  square  error. 
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II,  Statement  of  the  Pro,  'em 

| 

Throughout  the  paper,  we  assume  for  convenience  that  we  are 
dealing  with  incoherent  illumination.  All  observed  and  estimated 
quantities  are  intensities,  and  the  transfer  function  of  the  optical 
system  is  that  appropriate  for  incoherent  light  [2],  This  assumption 
simplifies  the  discussion  somewhat  and  avoids  confusion,  but  it  is  by 
no  means  essential. 

The  object  to  be  estimated  is  denoted  by  x{a,  p),  where  a  and  p 
are  the  object-plane  coordinates.  For  simplicity  of  presentation,  it 
is  assumed  that  the  optical  system  exhibits  spatial  invariance  [  3],  and 
that  it  has  a  point  spread  function  a{|,  tj)  where  £  and  q  are  the  image- 
plane  coordinates.  The  configuration  is  illustrated  in  Figure  1, 

The  image  of  x(a,  p)  is  given  by 

f 

y(i.  n>  «  'll*  a(5,  r,)  (l) 

1 

where  x(£,  tj)  is  the  object  distribution  referred  to  the  image  plane  and 

I 

the  star  denotes  convolution.  This  image  is  corrupted  by  additive  noise 

| 

n{£,  h)*  an<*  t*16  resultant  image  z(|,  r|)  =  y(£,  tj)  +  n(|,  q)  is  available 
for  processing.  The  objective  is  to  estimate  x(a,  p)  based  on  the  ob- 


gorvAiion  c(£,  t})  and  any  a  priori  information  available  about  the 
object. 


HI.  Sampling  and  Vector  Representation 

Almost  any  reasonable  object  x(a(  p)  can  be  represented  to  the 
desired  degree  of  accuracy  by  partitioning  the  object  plane  into  small 
regions  over  which  x(a,  fi)  is  essentially  constant.  If  these  regions  are 
made  sufficiently  small,  each  of  them  can  be  treated  as  an  impulse 
or  point  source  of  light.  The  image  due  to  one  of  these  regions  is 
then  approximately  proportional  to  the  point  spread  function  of  the 
optical  system  displaced  according  to  the  location  of  that  region.  The 
total  image  is  the  superposition  of  ihe  images  of  all  the  elementary 
regions. 

Let  the  intensity  of  the  k«th  region  in  the  object  plane  be  s 

X^ak*  ^k^'  i^^g®  °*  this  region  is  approximately  a(£  -  4^, 

ti  -  tv  ),  where  £.  and  tu  are  the  image-plane  coordinates  corresponding 
k  k  k 

to  c,j^  and  (3^.  If  the  object  consists  of  n  such  regions,  the  total  image 


is 


n 


y(S,  t})  =  L  a(|  -  ,  r\  -  tik) 


(2) 


i 
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If  thin  image  is  sampli  i  at  the  points  (|  , 


n  )»••«» 


m 

h 


>, 


we  have 


y(4J»  -i) 


n 


k=l 


k* 


J  v 

‘1  ~v 


(3) 


Note  that  t^)  is  not  necessarily  equal  to  (c  r(.). 

J  J 

a^k  »  a(4J  -  i-k,  r,J  -  qk)  and  =  y(£J,  qJ),  we  have 


If  we  set 


Yi 


k=l 


xi  a 

k  jk 


(4) 


It  matrix  form. 


(5) 


where 


The  samples  y.  are  in  gene  re  1  perturbed  by  noise  uampies  n  ,  and 
J  J 


HmwmiWiwiiH 
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the  observed  vector  is 


=  jr+n  =  Ax+n 


where  n  =  /  .  1 


Throughout  the  remainder  of  the  paper,  it  is  assumed  that  the  approximate 
vector  representation  discussed  above  is  sufficiently  accurate  that  the 
associated  error  can  be  neglected.  Heme,  the  vector  x  will  be  regarded 
as  the  object,  and  the  estimation  procedures  considered  will  be  evaluated 
in  terms  of  how  closely  they  reproduce  x. 


IV.  The  Optimum  Estimation  Procedure 


In  order  to  proceed  further,  we  now  make  some  additional  assumptions 

about  the  noise  and  the  a  priori  infe  rmatiorv.  We  assume  that  the  noise 

vector  n  has  zero  mean  anu  known  covariance  matrix  K  .  The  restriction 

— n 

to  zero  mean  is  unimportant,  siir  e  only  the  fl actuations  about  the  mean 
limit  system  performance.  Further,  we  assume  that  the  a  priori 
information  about  the  object  to  be  estimated  consists  of  the  mean  vector 
£  and  the  covariance  matrix  and  that  x  and  n  are  independent.  The 
quantities  g  and  do  not  necessarily  correspond  to  any  actual  statistical 
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fluctuations  which  the  object  undergoes,  but  rather  reflect  our  initial 
ignorance  of  the  object.  Roughly  speaking,  a  covariance  matrix  K  with 

— X 

small  entries  corresponds  to  a  situation  in  which  we  have  considerable 
prior  infcrmation  about  the  object,  while  a  covariance  matrix  with 
large  entries  corresponds  to  a  situation  in  which  we  have  little  prior 
information. 

The  approach  to  the  estimation  problem  taken  in  this  paper  is  that 

of  discrete  Wiener  filtering  of  the  perturbed  image  samples.  That  is, 

we  seek  that  discrete  linear  filter  whose  output  x  minimizes  the  mean 

a  A 

square  error  E[(x  -  x)'  ^x  -  x)] ,  where  the  prime  denotes  vector  or 
matrix  transpose. 

We  first  assume  that  the  prior  meanjz  is  equal  to  zero.  The  linear 
filtering  operation  on  the  observed  vector  z  is  then  represented  by  matrix 
multiplication,  and  any  linear  estimate  of  x  is  a  vector  of  the  form 
Hz.  Our  objective  is  to  choose  the  matrix  H  to  minimize  the  mean 
square  error 


3  “  E  *  -  x)1  (H  z  -  x)  ]  (7) 

Using  the  fact  that  z  =  Ax  +  n,  this  expression  can  be  written  in  the 


form 
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©  =  E  {[  (HA  -  I_)  x]  ’  [  (HA-I)x] 

+  [  (IJA  -I)x]'Hn+  {Hn )•  (HA  -  ijx 
+  (Hn)1  Hji  ^ 

Since  x  and  n  are  assumed  to  be  independent  random  vectors  with 
zero  means,  the  cross-product  terms  in  (8)  are  zero  and  the  error 
expression  reduces  to 

e  =  E[[  (HA  -  J_)  x]«  [  (HA  -  _I)  x]  +  (H  n)»  Hn  j 
=  tr(  (HA-I)  Kx(HA-I)»  +  HKnH']  (9) 

where  tr  denotes  the  trace  of  a  matrix  [4j. 

Upon  expanding  (9)  and  completing  the  square,  we  obtain  an 
expression  for  e  which  can  be  minimized  by  inspection. 

I  i 

e  =  tr  f[H(AK  A«  +  K  )2  -  K  A'  (AK  A1  +  K  )“*] 

i  .1 

[  H(A K  A'  +  K  )2  -  K  A*  (AK  A«  +  K  }  *}* 

—  —  — x  —  — n  — x  —  —  — x  —  — n 

+  K  -  K  A«  (AK  A'  +  K  )_1  AK  ?  (10) 

—x  —x— - x  —  — n - x  J 

X 

where  (B)2  is  the  symmetric  square  root  of  the  symmetric  matrix  B 
[  5].  Since  the  term  involving  H  in  (10)  is  non-negative,  the  mean 
square  error  e  achieves  its  minimum  for  that  H  which  satisfies 
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i  1 

[  H (A  K  A*  +  K  )2  -  K  A'  (A  K  A*  +  K  )'*]  =  0 
— x—  — n  — x— x—  — n 


This  yields  as  the  optimum  discrete  linear  filter  for  the  zero-mean 
case 


(11) 


A  somewhat  more  convenient  form  for  our  purposes,  obtained  after 
soma  manipulation  of  matrices,  is 

H  -  (A«  K  ‘  1  A  +  K  '  V  1  A'  K  ‘  1  (12) 

—  —  — n  —  -x  — n 

Thus  far  we  have  assumed  that  the  prior  meanji  is  zero.  The 
above  results  are  now  extended  to  the  general  case  in  which  fx  4  0, 

This  extension  is  accomplished  by  noting  that  the  above  results  are 
applicable  to  estimating  the  amount  by  which  x  deviates  fromjx,  namely 
x  -  fi.  Hence,  we  merely  apply  the  above  result  to  x  -  jx  and  add£  to 
the  resulting  estimate  of  the  deviation.  This  yields  as  the  optimum 
linear  estimate  of  th  object  distribution 


A 

x  =  £  +  H(z  -  Ap) 

=  Hz  +  (I  -  HAfcx  (  a 3 ) 
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which  becomes  upon  substituting  (12)  for  H 

A  _  1  - 1  - 1  - 1 

x  =  (A'K  A  +  K  )  A»K  z 
—  —  — n  —  —x  -  —  n  — 

4  {A’K  A  +  K  _1f  1  K  ‘i  u 
—  — n  —  — -X  — 


(14) 


This  is  the  desired  general  result. 

V.  Special  Cases 

A 

It  is  seen  that  the  optimum  estimate  x  is  just  a  linearly- filtered 
combination  of  the  observed  vector  z  and  the  a  priori  meanjv,  The 
details  of  this  filtering  and  combining  depend  on  the  noise,  the  optical 
system,  and  the  a  priori  information.  Since  the  operations  involved 
in  obtaining  the  optimum  estimate  are  in  general  rsather  complicated, 
we  now  consider  several  special  cases  in  order  to  gain  additional 
insight  into  the  estimation  procedure, 

A.  K  =  <r  2  I,  K  =  «r  2  I 
— n  n  —  —x  x  — 

2  2 

To  begin  with,  we  assume  that  K  «=  <r  I  and  X  =  ir  I, 

—n  n  -  -x  x  — 

where  1^  is  the  identity  matrix.  This  corresponds  to  a  situation  in 
which  the  noise  is  white  and  in  which  information  about  one  coordinate 
of  x  provides  no  information  about  any  other  coordinate.  In  this  case,  the 


. . . . m . . . n . •mmm*'  * . . . 
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optimum  estimate  becomes 

x  * 

+ 

This  expression  provides  some  insight  regarding  the  relative 

weighting  of  the  observation  £  and  the  a  priori  mean  If  the  a  priori 

uncertainty  about  x  is  much  larger  than  the  uncertainty  due  to  the 

2  2 

additive  noise  (the  usual  case),  then  <r  >>  cr  .  From  (15)  we  see 

x  n 

that  in  this  case  the  optimum  estimation  procedure  weights  the  cb- 

serva.ion  z  much  more  heavily  than  the  a  priori  mean  p  provided  that 

the  elements  of  A  are  not  too  small.  This  agrees  with  intuition,  since 

if  the  a  priori  uncertainty  is  very  large  v/e  have  little  confidence  in  the 

2  2 

a  priori  mean.  On  the  other  hand,  if  cr^  >>  ,  we  have  a  very  noisy 

measurement.  From  (15),  we  see  that  in  this  case  j£  is  weighted  more 
heavily  than  z,  again  in  agreement  with  intuition,  and  the  experiment 
is  not  a  very  informative  one. 


A*  A  1  .-1  A» 

+  rr 


n 


n 


A*  A  1,-1  1 

f =-f=-  +  — -h— }  - ?  H 

n  x  x 


(15) 
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B.  Large  a  priori  uncertainty 

We  can  obtain  additional  insight  into  the  optimum  estimation 
procedure  by  looking  at  the  important  limiting  case  of  large  a  priori 
uncertainty.  The  optimum  estimate  in  this  case  is  obtained  by  making 
a  Neumann  series  expansion  [  6]  of  (14)  and  retaining  only  the  first 
texm.  This  procedure  yields 

x  *  <A‘K  “l  Af  1  A»K  _1  z  (16) 

—  —  —a  —  —  —  n  — 


This  procedure  is  now  considered  under  various  assumptions  about 

A  and  K  . 

—  — n 

If  the  system  matrix  A  is  square  and  invertible,  (16)  becomes 

x  =  (A*  K  "  1  A)~  1  A*  Kn"  1  AA~  1  z  =A'  1  z  (17) 


In  this  case,  the  optimum  estimation  procedure  merely  "divides  out" 

the  effects  of  the  optical  system  by  multiplying  the  observed  vector  by 
-  1 

A  .  This  is  precisely  the  procedure  which  would  be  followed  in  the 
absence  of  noise,  and  the  above  analysis  clearly  indicates  the  limitations 
on  this  procedure  and  the  conditions  under  which  it  is  optimum.  These 
conditions  are  negligible  a  priori  information  about  the  object  to  be 
estimated  (except  for  information  regarding  spatial  extent)  and  a 


. mmmmmrnm . mm  *  <  . mmmmmmmmim'  otmm! 
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sampling  scheme  in  the  image  plane  which  yields  a  square  and  in¬ 
vertible  system  matrbc  A.  If  either  of  these  conditions  is  violated, 
the  more  general  data  processing  procedures  given  by  (14)  or  (16) 
should  be  used* 

A  further  simplification  of  (16)  is  achieved  in  the  event  that  the 

additive  noise  is  white.  In  this  case,  K  =  <r  ^  I  and 

— n  n  — 


x  =  (A*  A)"1  JJ  z  (18) 

This  expression  reduces  to  (171  of  course,  if  A  is  square  and  invertible. 

VI.  Evaluation  of  the  Optimum  Procedure 

For  the  purpose  of  evaluating  the  optimum  estimation  or  restoration 

procedure,  we  assume  that  the  actual  (but  unknown)  object  vector  is 

x  .  The  observed  vector  is  thus 
— o 

z  =  A  x^  +  n 

We  further  assume  that  the  a  priori  uncertainty  is  very  large,  which 
leads  us  to  use  the  estimate  given  by  (16).  This  yields 
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x  =  (A1  K^~  1  A)”  1  A1  K^"  1  (A  xq  +  n) 

~  x  4  (A*  K  ”  1  A)”  1  A*  K  ” 1  n 
— o  -  -n  — - -n  - 


[19) 


The  mean  square  error  (conditioned  on  the  actual  object  vector  xq)  is 
given  by 


e  =  E[  (x  -  Xo)«  (x  -  xq)] 

=  E  f  [(^K/^r^'K/1  n]»  [(A'Kn"1A)"lAfK|‘1ri]| 


which  after  some  manipulation  of  matrices  can  be  written 


e  =  tr  [(A'Kn_1  Af1] 


In  the  special  case  of  white  noise.  K  «  <r  I  and 

— n  n  — 


e  -  o-^2  tr  [  (A*  A)'  l] 


(20) 


(21) 


Recall  that  the  terms  in  the  matrix  A  are  determined  by  the  point 
spread  function  of  the  optical  system  and  the  campling  points.  Similarly, 
the  covariance  matrix  is  determined  by  the  sampling  scheme  in  the 
image  plane  as  well  as  by  the  nature  of  the  noise.  Hence,  it  is  clear 
that  we  can  control  the  form  of  the  matrix  A*  K  *  A  io  a  limited  extent 


L 
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through  our  choice  of  sampling  schemes.  A.  judicious  choice  of  sampling 
points  can  reduce  the  mean  square  error  associated  with  the  optimum 
estimation  procedure. 

The  problem  of  where  to  sample  and  how  many  samples  to  take 
is  both  important  and  nontrivial,  and  is  currently  being  investigated 
Although  no  general  solution  is  available  at  this  point,  some  of  the 
features  of  this  solution  aro  apprent.  First,  it  is  clear  that  samples 
near  the  central  part  of  the  image  are  more  useful  than  samples  on 
the  fringes  where  the  signal- to- noise  ratio  is  low.  Second,  it  is  not 
difficult  to  show  that  adding  additional  samples  without  changing  the 
location  of  previous  samples  will  never  increase  the  error,  and  in  fact 
will  decrease  it  except  in  special  cases.  Thus,  it  may  well  be  that  we 
will  want  to  take  more  samples  than  we  have  unknowns  to  estimate 
(m  >  n).  The  usefulness  of  additional  samples  is  limited  by  the  increase 
in  the  complexity  of  the  estimation  procedure,  of  course,  since  the 
mum  estimate  for  th  "  case  is  given  by  (16)  where  A  is  an  m  x  n 
matrix  rather  than  by  the  intuitive  expression  (17)  with  A  an  n  x  n  matrix. 


\ 
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YU,  Example 


We  now  consider  a  simple  example  to  illustrate  several  of  the 
above  ideas.  The  configuration  for  this  example  is  shown  in  Figure  2. 
The  object  is  a  single  line  source  at  the  origin  in  the  object  plane,  and 
it  is  desired  to  estimate  the  intensity  of  this  source.  We  assume  in  all 
cases  that  the  a  priori  uncertainty  is  very  large.  The  aperture  is  a 
rectangular  slit  with  point  spread  function 

a(£)  =  a  sine2  £ 

2  , 

=  a  81-n-  (22) 

<nr 

The  noisy  image  is  then  given  by 

y(£)  =  xa  sinc2£  +  n{£)  (23) 

To  begin  with,  we  consider  the  case  where  a  single  sample  is 
taken  at  £  =  0.  We  then  have 


z  = 


ax  +  n 


(24) 


The  matrix  A  in  this  case  rednres  to  the  scalar  a,  and  the  optimum 
estimate  is,  from  (17) 


A 

x  =  z/a 


(25) 


The  mean  square  error  associated  with  this  estimate  is,  from  (21), 


e  = 


n 


(26) 


2  . 

where  <r^  is  the  variance  of  the  noise  sample  at  £  =0, 

Now  suppose  that  in  addition  to  the  sample  at  £  =  0,  a  second 

sample  is  taken  at  £  =  1/2,  Assume  for  the  moment  that  K  =  tr  ^  I; 

— n  n  — 

i*e. ,  the  noise  samples  are  uncorrelated  and  have  the  same  variance. 
The  system  matrix  A  in  this  case  becomes 


a 


-  =  ^a/ir^  =  *0.404  a* 


Since  A  is  no  longer  square,  (16)  rather  than  (17)  must  be  used  to 
estimate  x.  The  optimum  estimate  for  this  case  is 


x  =  0.  859  Zj/a  +  0.  347z^,/a 


(27) 


where  is  the  observed  sample  at  the  origin  and  z^that  at  £  =  1/2. 
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The  mean  square  error  associated  with  this  estimate  is  seen  from 
(21)  to  be 

e  =  0.  859  r  ^/a^,  (28) 


a  moderate  reduction  over  that  associated  with  the  single- sample 
procedure. 

To  illustrate  the  effects  of  correlation  between  noise  samples 
in  a  particular  case,  we  consider  the  same  two  sampling  points  as 
above,  but  now  assume 


2 

<r 

n 


The  optimum  estimate  is  seen  from  (16)  to  be 


x  *  1. 052  z  /a  .  0.  127  z2/a 


(29) 


with  an  associated  mean  square  error 


e 


0.  989 


2,  2 

<r  /a 


n 


(30) 


The  effects  of  correlation  between  noise  samples  in  this  example  may¬ 
be  assessed  by  comparing  (29)  and  (30)  with  (27)  and  (28). 
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Vm.  Multiple  Olservations 


In  many  situations,  we  make  not  one  but  several  observations 
of  the  object  to  be  estimated*  For  example,  we  may  have  several 
different  photographs  of  the  object.  If  the  noise  samples  associated 
with  different  observations  a^'e  independent,  we  can  make  use  of  this 
fact  to  improve  our  estimate. 

1  Q 

Suppose  that  q  observations  z  .  .  z  are  made,  where 

i  »  i 

z  =  “  x  +  n 

The  n1  are  assumed  to  independent  noise  vectors,  each  having 
zero  mean  and  covariance  matrix  K^.  The  optimum  estimate  based 
on  the  observations  7}  .  ,  z^  can  be  shown  by  a  method  similar  to  that 
used  for  a  single  observation  to  be 

x  =  (qA»K  _1  A  +  K  ~l)~l  qA'K  “lz 
— q  —  — n  —  —x  -  — n  — 

+  (iA«K“1A  +  K"lr1K“1tt  (31) 

—  — n  —  —x  —x  — 

where 

1  q  i 

z  =  -  S  z  (32) 

"  *  i=l” 
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is  the  sample  mean  of  the  observed  vectors.  This  estimate  is  identical 
in  form  to  that  given  by  (14)  for  a  single  observation,  and  it  reduces 
to  (14)  for  q  =  1  as  it  should* 

If  the,  matrix  A*  *  A  is  non- singular,  as  it  should  be  for  good 
restoration,  we  can  consider  the  limiting  value  of  the  estimate  as 
the  number  of  observations  becomes  arbitrarily  large.  We  have 


lim 
q  -*■  oo 


A 

X 

-q 


(A1  K  -lA)~lA*K  '  1  z 
—  — n  —  —  — n  — 


(33) 


But  by  the  law  of  large  numbers  [7],  we  further  have 


lim  _ 
z 

q  -►ao  — 


Ax 
- — 0 


(34) 


with  probability  one,  where  x^  is  the  true  value  of  the  object,  assumed 
to  be  constant  throughout  the  experiment.  The  optimum  estimate  aus 
reduces  in  the  limit  to 


lim  a 
x 

q  —  oo  — q 


(A‘K  "1A)“lA«K  _1Ax  =  x  (35) 
--  —  n  —  —  — n  — o  —o 


and  we  see  that  the  object  is  reproduced  exactly  despite  the  presence 
of  noise. 


1 

Object 


Aperture 


Figure  2.  Lrr.age  of  a  line  source  through  a  rectangular  aperture. 
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